home *** CD-ROM | disk | FTP | other *** search
/ The Very Best of Atari Inside / The Very Best of Atari Inside 1.iso / mint / mntlib25 / _muldf3.cpp < prev    next >
Text File  |  1992-12-12  |  4KB  |  189 lines
  1. | mjr: not needed on the TT
  2.  
  3. #ifndef    __M68881__
  4.     .text
  5.     .even
  6.     .globl    __muldf3, ___muldf3
  7.  
  8. __muldf3:
  9. ___muldf3:
  10.  
  11. # ifdef    sfp004
  12.  
  13. | double precision floating point stuff for Atari-gcc using the SFP004
  14. | developed with gas
  15. |
  16. | double precision multiplication
  17. |
  18. | M. Ritzert (mjr at dfg.dbp.de)
  19. |
  20. | 4.10.1990
  21. |
  22. | no NAN checking implemented since the 68881 treats this situation "correctly",
  23. | i.e. according to IEEE
  24.  
  25. | addresses of the 68881 data port. This choice is fastest when much data is
  26. | transferred between the two processors.
  27.  
  28. comm =     -6
  29. resp =    -16
  30. zahl =      0
  31.  
  32. | waiting loop ...
  33. |
  34. | wait:
  35. | ww:    cmpiw    #0x8900,a0@(resp)
  36. |     beq    ww
  37. | is coded directly by
  38. |    .long    0x0c688900, 0xfff067f8
  39.  
  40.     lea    0xfffffa50:w,a0
  41.     movew    #0x5400,a0@(comm)    | load first argument to fp0
  42.     cmpiw    #0x8900,a0@(resp)    | check
  43.     movel    a7@(4),a0@
  44.     movel    a7@(8),a0@
  45.     movew    #0x5423,a0@(comm)
  46.     .long    0x0c688900, 0xfff067f8
  47.     movel    a7@(12),a0@
  48.     movel    a7@(16),a0@
  49.     movew    #0x7400,a0@(comm)    | result to d0/d1
  50.     .long    0x0c688900, 0xfff067f8
  51.     movel    a0@,d0
  52.     movel    a0@,d1
  53.     rts
  54.  
  55. # else    sfp004
  56.  
  57. | double floating point multiplication routine
  58. |
  59. | written by Kai-Uwe Bloem (I5110401@dbstu1.bitnet).
  60. | Based on a 80x86 floating point packet from comp.os.minix, written by P.Housel
  61. |
  62. |
  63. | Revision 1.2, kub 01-90 :
  64. | added support for denormalized numbers
  65. |
  66. | Revision 1.1, kub 12-89 :
  67. | Ported over to 68k assembler
  68. |
  69. | Revision 1.0:
  70. | original 8088 code from P.S.Housel
  71.  
  72. BIAS8    =    0x3FF-1
  73.  
  74.     lea    sp@(4),a0
  75.     moveml    d2-d7,sp@-
  76.     moveml    a0@,d4-d5/d6-d7 | d4-d5 = v, d6-d7 = u
  77.     subw    #16,sp        | multiplication accumulator
  78.  
  79.     movel    d6,d0        | d0 = u.exp
  80.     swap    d0
  81.     movew    d0,d2        | d2 = u.sign
  82.     lsrw    #4,d0
  83.     andw    #0x07ff,d0    | kill sign bit
  84.  
  85.     movel    d4,d1        | d1 = v.exp
  86.     swap    d1
  87.     eorw    d1,d2        | d2 = u.sign ^ v.sign (in bit 31)
  88.     lsrw    #4,d1
  89.     andw    #0x07ff,d1    | kill sign bit
  90.  
  91.     andl    #0x0fffff,d6    | remove exponent from u.mantissa
  92.     tstw    d0        | check for zero exponent - no leading "1"
  93.     beq    0f
  94.     orl    #0x100000,d6    | restore implied leading "1"
  95.     bra    1f
  96. 0:    addw    #1,d0        | "normalize" exponent
  97. 1:    movel    d6,d3
  98.     orl    d7,d3
  99.     beq    retz        | multiplying by zero
  100.  
  101.     andl    #0x0fffff,d4    | remove exponent from v.mantissa
  102.     tstw    d1        | check for zero exponent - no leading "1"
  103.     beq    0f
  104.     orl    #0x100000,d4    | restore implied leading "1"
  105.     bra    1f
  106. 0:    addw    #1,d1        | "normalize" exponent
  107. 1:    movel    d4,d3
  108.     orl    d5,d3
  109.     beq    retz        | multiplying by zero
  110.  
  111.     addw    d1,d0        | add exponents,
  112.     subw    #BIAS8+16-11,d0    | remove excess bias, acnt for repositioning
  113.  
  114.     clrl    sp@        | initialize 128-bit product to zero
  115.     clrl    sp@(4)
  116.     clrl    sp@(8)
  117.     clrl    sp@(12)
  118.     lea    sp@(8),a1    | accumulator pointer
  119.  
  120. | see Knuth, Seminumerical Algorithms, section 4.3. algorithm M
  121.  
  122.     swap    d2
  123.     movew    #4-1,d2
  124. 1:
  125.     movew    d5,d3
  126.     mulu    d7,d3        | mulitply with bigit from multiplier
  127.     addl    d3,a1@(4)    | store into result
  128.     movew    d4,d3
  129.     mulu    d7,d3
  130.     movel    a1@,d1        | add to result
  131.     addxl    d3,d1
  132.     movel    d1,a1@
  133.     roxlw    a1@-        | rotate carry in
  134.  
  135.     movel    d5,d3
  136.     swap    d3
  137.     mulu    d7,d3
  138.     addl    d3,a1@(4)    | add to result
  139.     movel    d4,d3
  140.     swap    d3
  141.     mulu    d7,d3
  142.     movel    a1@,d1        | add to result
  143.     addxl    d3,d1
  144.     movel    d1,a1@
  145.  
  146.     movew    d6,d7
  147.     swap    d6
  148.     swap    d7
  149.     dbra    d2,1b
  150.  
  151.     swap    d2        | [TOP 16 BITS SHOULD BE ZERO !]
  152.  
  153.     moveml    sp@(2),d4-d7    | get the 112 valid bits
  154.     clrw    d7        | (pad to 128)
  155. 2:
  156.     cmpl    #0x0000ffff,d4    | multiply (shift) until
  157.     bhi    3f        |  1 in upper 16 result bits
  158.     cmpw    #9,d0        | give up for denormalized numbers
  159.     ble    3f
  160.     swap    d4        | (we''re getting here only when multiplying
  161.     swap    d5        |  with a denormalized number; there''s an
  162.     swap    d6        |  eventual loss of 4 bits in the rounding
  163.     swap    d7        |  byte -- what a pity 8-)
  164.     movew    d5,d4
  165.     movew    d6,d5
  166.     movew    d7,d6
  167.     clrw    d7
  168.     subw    #16,d0        | decrement exponent
  169.     bra    2b
  170. 3:
  171.     movel    d6,d1        | get rounding bits
  172.     roll    #8,d1
  173.     movel    d1,d3        | see if sticky bit should be set
  174.     orl    d7,d3        | (lower 16 bits of d7 are guaranteed to be 0)
  175.     andl    #0xffffff00,d3
  176.     beq    4f
  177.     orb    #1,d1        | set "sticky bit" if any low-order set
  178. 4:    addw    #16,sp        | remove accumulator from stack
  179.     jmp    norm_df        | (result in d4/d5)
  180.  
  181. retz:    clrl    d0        | save zero as result
  182.     clrl    d1
  183.     addw    #16,sp
  184.     moveml    sp@+,d2-d7
  185.     rts            | no normalizing neccessary
  186.  
  187. # endif    sfp004
  188. #endif    __M68881__
  189.